/* ----------------------------------------------------------------- */
/*  Copyright (c) 2014                                               */
/*  Fernando Noveletto Candiani         7239131                      */
/*                                                                   */
/*  This program is free software; you can redistribute it and/or    */
/*  modify it under the terms of the GNU General Public License as   */
/*  published by the Free Software Foundation; either version 3 of   */
/*  the License, or (at your option) any later version.  See the     */
/*  file LICENSE included with this distribution for more            */
/*  information.                                                     */
/*  email: fncandiani@usp.br                                         */
/* ----------------------------------------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>

#define _AIR_DENSITY_ 1.2922 //In Kg/m^3
#define _AIR_CONSTANT_ -0.000101741 //Result of the hypsometric equation
#define _EARTH_RADIUS_ 6335437.0 //In meters
#define _SEE_LEVEL_GRAVITY_ 9.80665 // In m/s^2

#define _ANGLE_INCREMENT_ 1E-02 // The specified error is 10^-3

#define _RADIUS_INCREMENT_ 1E-03 //
#define _RADIUS_MAX_ 10E06

#define _H_ 1E-04 //Integration step

typedef struct Params Params;

struct Params{
	double m_0, m_p, t_q, F_p, alpha, Cd, r; //constants read from stdin (r is simulated as being read from stdin)
	double height, reach, angle, speed, time; // variables that changes throughout the time
	double x, y; // variables that will be used in the runge kutta method
	int offset;
	double result;
};

pthread_t *callThd;
int num_threads;
char option;
double S_f;

//Calculates the gravity that is applied to the rocket (m/s^2)
double gravity(double height){
	return (_SEE_LEVEL_GRAVITY_)*((_EARTH_RADIUS_*_EARTH_RADIUS_)/((_EARTH_RADIUS_+height)*(_EARTH_RADIUS_+height)));
}

//Calculates the drag force that is applied to the rocket
double drag_force(double aerodynamic_drag, double cross_section, double air_density, double speed){
	return (aerodynamic_drag*cross_section*air_density*speed*speed)/2.0f;
}

//Calculates the air density and returns it Kg/m^3
double air_density(double height){
	return _AIR_DENSITY_*exp(_AIR_CONSTANT_*height);
}

//Calculates the mass and its buoyancy of the rocket in the time t
void rocket_mass(double total_mass, double fuel_mass, double time, double fuel_time, double buoyancy_fuel, double *mass, double *buoyancy){

	if(time > fuel_time){
		*mass = total_mass - fuel_mass;
		*buoyancy = 0.0f;
	}else{
		*mass = total_mass - (fuel_mass*(time/fuel_time));
		*buoyancy = buoyancy_fuel;
	}

}

//Returns the rocket cross section assuming it is a cylinder
double rocket_cross_section(double radius){
	return M_PI * radius * radius;
}

//dV/dt = ...
double speed_func(Params par){
	double rkt_mass, buoyancy;

	//Calculates the mass and the buoyancy force
	rocket_mass(par.m_0, par.m_p, par.y, par.t_q, par.F_p, &rkt_mass, &buoyancy);

	//Calculates the speed function defined in the description

	return ((buoyancy*cos(par.alpha))/rkt_mass) -
			(drag_force(par.Cd, rocket_cross_section(par.r), air_density(par.height), par.x)/rkt_mass) -
			(gravity(par.height)*sin(par.angle));
}

//d gamma/dt = ...
double angle_func(Params par){

	if(0.0f == par.speed){
		return 0.0f;
	}

	double rkt_mass, buoyancy;

	//Calculates the mass and the buoyancy force
	rocket_mass(par.m_0, par.m_p, par.y, par.t_q, par.F_p, &rkt_mass, &buoyancy);

	return ((buoyancy*sin(par.alpha))/(rkt_mass*par.speed)) +
		   ( gravity(par.height)*cos(par.x)/par.speed);
}

//dS/dt = ..
double reach_func(Params par){
	return par.speed*cos(par.angle);
}

//dH/dt = ..
double height_func(Params par){
	return par.speed*sin(par.angle);
}

double calculate_k1(Params par, double(*function)(Params)){
	return function(par);
}

double calculate_k2(double h, double k1, Params par, double(*function)(Params)){
	par.x = par.x + (h/2);
	par.y = par.y + ((h/2)*k1);
	return function(par);
}

double calculate_k3(double h, double k2, Params par, double(*function)(Params)){
	par.x = par.x + (h/2);
	par.y = par.y + ((h/2)*k2);
	return function(par);
}

double calculate_k4(double h, double k3, Params par, double(*function)(Params)){
	par.x = par.x + h;
	par.y = par.y + h*k3;
	return function(par);
}

double rungekutta4(double h, Params par,  double(*function)(Params)){
	double k1 = calculate_k1(par,function);
	double k2 = calculate_k2(h, k1, par, function);
	double k3 = calculate_k3(h, k2, par, function);
	double k4 = calculate_k4(h, k3, par, function);
/*	printf("k1 %lf \tk2 %lf \tk3 %lf \tk4 %lf \tr: %lf\n", k1,k2,k3,k4, par.x + ( h/6 )*( k1+2*(k2+k3)+k4));*/
	return par.x + ( h/6 )*( k1+2*(k2+k3)+k4);
}


double rungekutta4_wrapper(Params par){

	//The rocket reached the ground and it is not the start
	while(!(0.0f >= par.height && 0.1f <= par.time)){

		//Calculates the new speed
		par.x = par.speed;
		par.y = par.time;
/*		printf("speed(V) -- \t");*/
		double speed = rungekutta4(_H_, par, speed_func);

		//Calculates the new angle
		par.x = par.angle;
		par.y = par.time;
/*		printf("angle(Y) -- \t");*/
		double angle = rungekutta4(_H_, par, angle_func);

		//Calculates the new reach
		par.x = par.reach;
		par.y = par.time;
/*		printf("reach(S) -- \t");*/
		double reach = rungekutta4(_H_, par, reach_func);

		//Calculates the new height
		par.x = par.height;
		par.y = par.time;
/*		printf("height(H) -- \t");*/
		double height = rungekutta4(_H_, par, height_func);

		par.time = par.time + _H_;
		par.speed = speed;
		par.angle = angle;
		par.reach = reach;
		par.height = height;
	}
	return par.reach;
}

void* rungekutta4_thd(void *arg){
	Params *par = (Params *) arg;
	double angle = 0.0f;
	double first_angle = par->offset*(M_PI/num_threads);
	double final_angle = (par->offset+1)*(M_PI/num_threads) - _ANGLE_INCREMENT_;

	if('G' == option){
		for(angle = first_angle; angle < final_angle; angle += _ANGLE_INCREMENT_){
			par->angle = angle;
/*			par->angle = 0.69813f;*/
			par->speed = 0.0f;
			par->time = 0.0f;
			par->height = 0.0f;
			par->reach = 0.0f;
			double S = rungekutta4_wrapper(*par);
			if(S == S_f){
				par->result = angle;
				break;
			}
		}
	}else{
		for(par->r = _RADIUS_INCREMENT_; par->r < _RADIUS_MAX_; par->r += _RADIUS_INCREMENT_){
			for(angle = first_angle; angle < final_angle; angle += _ANGLE_INCREMENT_){
				par->angle = angle;
				par->speed = 0.0f;
				par->time = 0.0f;
				par->height = 0.0f;
				par->reach = 0.0f;
				double S = rungekutta4_wrapper(*par);
				if(S == S_f)
					par->result = angle + par->offset;
			}
		}
	}

	pthread_exit((void*) EXIT_SUCCESS);
}

void cpy_params(Params par, Params *n_par){
		(*n_par).m_0 = par.m_0;
		(*n_par).m_p = par.m_p;
		(*n_par).t_q = par.t_q;
		(*n_par).F_p = par.F_p;
		(*n_par).alpha = par.alpha;
		(*n_par).Cd = par.Cd;
		(*n_par).r = par.r;
}


int main(int argc, char **argv){

	if(argc < 2){
		num_threads = 4;
	}else{
		num_threads = atoi(argv[1]);
	}

	int i;
	Params par[num_threads+1];

	scanf("%s %lf %lf %lf %lf %lf %lf", &option, &par[0].m_0, &par[0].m_p, &par[0].t_q, &par[0].F_p, &par[0].alpha, &par[0].Cd);

	if('G' == option){
		scanf("%lf", &par[0].r);
	}else{//Generates a NaN
		par[0].r = 0.0f/0.0f;
	}

	scanf("%lf", &S_f);

/*	printf("received - %c %lf %lf %lf %lf %lf %lf %lf %lf\n", option, par[0].m_0, par[0].m_p, par[0].t_q, par[0].F_p, par[0].alpha, par[0].Cd, par[0].r, S_f);*/

	callThd = (pthread_t *)malloc(sizeof(pthread_t)*num_threads);

	for (i = 0; i < num_threads; i++) {
		cpy_params(par[0], par+i+1);
		(par+i+1)->offset = i;
		(par+i+1)->result = 0.0f;
		pthread_create(&callThd[i], NULL, rungekutta4_thd, (void *) (par+i+1));
	}

	for (i = 0; i < num_threads; i++) {
		pthread_join(callThd[i], NULL);
		if( 0.0f == ((par+i+1)->result)){ // The result was not found by this thread
		}
	}

/*	printf("%c %.5lf %.5lf\n", option, gamma, r);*/
	printf("%.5lf %.5lf\n", (par+i)->result, par[0].r);

	return EXIT_SUCCESS;
}
